#8305211220 赵赟乾 
path="C:\\Users\\29228\\Desktop\\2024春学期各项汇报\\R语言与生物信息学\\8305211220 赵赟乾 R语言大作业"
setwd(path)



#数据分布直方图#####
library(readxl)  
# 使用read_excel函数读取Excel文件  
df = read_excel("GSE224239_genes_fpkm_expression.xlsx")  
library(tidyr)  
library(dplyr)
library(tidyverse)
library(pheatmap)
library(ggplot2)
#WT<-df[,-(3:21)]
DF<-df[,-1]
DF1<-DF[,-2]

write.csv(DF1,file="Data.csv")
# 读取csv文件
data <- read.csv("Data.csv")
data<-data[,-1]
# 检查重复的行名
duplicates <- data$gene_name[duplicated(data$gene_name)]
print(duplicates)  # 打印重复的行名

# 为每个重复的行名添加唯一标识符
data$unique_id <- make.unique(as.character(data$gene_name))

# 将新列 unique_id 设置为行名，并移除原来的 gene_name 列和 unique_id 列
rownames(data) <- data$unique_id
data <- data[,-which(names(data) == "gene_name")]
data<-data[,-ncol(data)]
# 检查结果
View(data)
write.csv(data,file="Data1.csv")
data<-as.matrix(data)
hist(data)
hist(data[data!=0])
hist(log2(data+1))
hist(log2(data[data!=0]))#取出0之后的数据对数化


#聚类分析#####
# 读取数据并进行对数转换

data1 <- read.csv("Data1.csv", row.names = 1)

data1 <- log2(data1 + 1) # 对数转换

# 数据转置
data2 <- t(data1) # 对样本进行聚类需要转置，把行设置为样本

# 计算样本间的欧式距离
dist <- dist(data2) # 欧式距离计算样本间距离

# 层次聚类
da <- hclust(dist, method = 'average') # 使用平均法进行层次聚类

# 打开 PDF 设备并绘制树状图
pdf('data1tree.pdf') # 设置保存的文件名
plot(da, hang = 0.5, cex = 0.4) # 绘制树状图，hang 控制标签位置，cex 控制标签大小
dev.off() # 关闭 PDF 设备




#差异分析######
#WT_CTRL：野生型小鼠的对照组，未接受LPS处理。
#WT_LPS：野生型小鼠经过LPS处理的组，用于观察LPS对野生型小鼠的炎症反应。
#KO_CTRL：Otud6a基因敲除小鼠的对照组，未接受LPS处理。
#KO-LPS：Otud6a基因敲除小鼠经过LPS处理的组，
#用count数据进行差异分析


# 加载必要的包
library(DESeq2)
# 读取csv文件，将第一列作为行名
data3 <- read.csv("Data1.csv", row.names = 1)
data4<-data3[,-(1:18)]
group <- substring(colnames(data4),10,11)#提取分组信息
group <- factor(group,labels = c("CTRL","LPS")) #分组信息改成CTRL和LPS,并变成因子型
coldata <- data.frame(row.names = colnames(data4), group) #建立meta信息
dds <- DESeqDataSetFromMatrix(countData = data4,
                              colData = coldata,
                              design = ~group)
dds$group<- relevel(dds$group, ref = "CTRL") # 指定CTRL作为对照组
dds2 <- DESeq(dds)
result <- results(dds2)
summary(result)
write.csv(result, file= "DESeq2结果.csv")

#火山图#####

# 加载必要的包
if (!requireNamespace("ggplot2", quietly = TRUE)) install.packages("ggplot2")
if (!requireNamespace("ggrepel", quietly = TRUE)) install.packages("ggrepel")
if (!requireNamespace("ggprism", quietly = TRUE)) install.packages("ggprism")
if (!requireNamespace("ggthemes", quietly = TRUE)) install.packages("ggthemes")

library(ggplot2)
library(ggrepel)
library(ggprism)
library(ggthemes)

# 读取CSV文件
plotdata <- read.csv("DESeq2结果.csv")
colnames(plotdata)[1] <- "ID"
view(plotdata)

# 修改列名
colnames(plotdata)[which(colnames(plotdata) == "log2FoldChange")] <- "log2fc"
colnames(plotdata)[which(colnames(plotdata) == "padj")] <- "P"

# 设置阈值
fccutoff <- 2  # fc阈值
pcutoff <- 0.05  # p阈值

# 新建一列group 全部设置为nosig
plotdata$group <- "nosig"

# 把满足条件的设置为up
plotdata$group[plotdata$log2fc > log2(fccutoff) & plotdata$P < pcutoff] <- "up"

# 把满足条件的设置为down
plotdata$group[plotdata$log2fc < -log2(fccutoff) & plotdata$P < pcutoff] <- "down"

# 新建一列label 用于图上标记信息
plotdata$label <- plotdata$ID
plotdata$label[plotdata$group == "nosig"] <- ""  # nosig不标记 设置为""

# 移除P值为0的行
plotdata <- plotdata[plotdata$P != 0, ]


# 绘制火山图函数
plot_volcano <- function(data, ylim = NULL, filename = NULL) {
  p <- ggplot(data, aes(x = log2fc, y = -log10(P))) +
    geom_point(aes(color = group)) +
    scale_color_manual(values = c("green", "grey", "red")) +
    geom_hline(yintercept = -log10(pcutoff), lty = 2) +  # 画横虚线
    geom_vline(xintercept = c(log2(fccutoff), -log2(fccutoff)), lty = 2) +  # 画竖虚线
    geom_text_repel(aes(label = label)) +
   
    theme_prism()
  
  if (!is.null(ylim)) {
    p <- p + scale_y_continuous(limits = ylim)
  }
  
  print(p)
  
  if (!is.null(filename)) {
    ggsave(filename, plot = p, width = 10, height = 12, units = "in", dpi = 300)
  }
}

# 显示火山图
plot_volcano(plotdata)
# 保存火山图为PNG文件并设置y轴限制
plot_volcano(plotdata, ylim = c(0, 350), filename = "v.png")



#富集分析#####
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

# 安装org.Mm.eg.db
BiocManager::install("org.Mm.eg.db")

# 安装DOSE, clusterProfiler和enrichplot
BiocManager::install(c("DOSE", "clusterProfiler", "enrichplot"))

# 安装ggplot2, ggrepel, ggprism和ggthemes
install.packages(c("ggplot2", "ggrepel", "ggprism", "ggthemes"))

# 加载包
library(org.Mm.eg.db)
library(DOSE)
library(clusterProfiler)
library(enrichplot)
library(ggplot2)
library(ggrepel)
library(ggprism)
library(ggthemes)

# 提取显著上调和下调的基因
up_genes <- plotdata$ID[plotdata$group == "up"]
down_genes <- plotdata$ID[plotdata$group == "down"]

# 将基因ID转换为ENTREZ ID
up_genes_entrez <- bitr(up_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Mm.eg.db)
down_genes_entrez <- bitr(down_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Mm.eg.db)
# GO分析
go_up <- enrichGO(gene = up_genes_entrez$ENTREZID,
                  OrgDb = org.Mm.eg.db,
                  ont = "BP", # BP: Biological Process, CC: Cellular Component, MF: Molecular Function
                  pAdjustMethod = "BH",
                  pvalueCutoff = 0.05,
                  qvalueCutoff = 0.2,
                  readable = TRUE)

go_down <- enrichGO(gene = down_genes_entrez$ENTREZID,
                    OrgDb = org.Mm.eg.db,
                    ont = "BP",
                    pAdjustMethod = "BH",
                    pvalueCutoff = 0.05,
                    qvalueCutoff = 0.2,
                    readable = TRUE)
# GO分析结果可视化
dotplot(go_up, showCategory = 10, title = "GO Enrichment Analysis (Upregulated Genes)")
ggsave("go_up_dotplot.png", plot = last_plot(), width = 10, height = 8, dpi = 300)
dotplot(go_down, showCategory = 10, title = "GO Enrichment Analysis (Downregulated Genes)")
# 保存结果

ggsave("go_down_dotplot.png", plot = last_plot(), width = 10, height = 8, dpi = 300)


#AVONA方差分析######
# 加载必要的R包
if (!requireNamespace("tidyverse", quietly = TRUE)) install.packages("tidyverse")
library(tidyverse)

# 读取数据
Data0 <- read.csv("Data1.csv", row.names = 1)

# 移除前12列
Data <- Data0[, -(1:12)]

# 将0替换为NA
Data[Data == 0] <- NA

# 对数变换
Data <- log2(Data)

# 提取样本名并分组
group <- colnames(Data)
group <- gsub("\\d", "", group)  # 去掉数字字符，即为分组信息

# 转置数据
testdata <- t(Data)

# 初始化一个空的p对象
p <- vector("numeric", ncol(testdata))

# 进行ANOVA方差分析
for (i in 1:ncol(testdata)) {
  value <- testdata[, i]
  
  # 提取各组不为NA的位置
  pos1 <- !is.na(value[group == "count.KO_CTRL"])
  pos2 <- !is.na(value[group == "count.KO_LPS"])
  pos3 <- !is.na(value[group == "count.WT_CTRL"])
  pos4 <- !is.na(value[group == "count.WT_LPS"])
  
  # 检查每组不为NA的数是否>=3
  if (sum(pos1) < 3 | sum(pos2) < 3 | sum(pos3) < 3 | sum(pos4) < 3) {
    p[i] <- NA
  } else {
    # 进行单因素方差分析
    aovtest <- oneway.test(value ~ group, var.equal = TRUE)
    p[i] <- aovtest$p.value
  }
}

# 创建数据框并保存结果
stat <- data.frame(ID = colnames(testdata), pvalue = p)

# 保存结果到文件
write.csv(stat, file = "anova_test_result.csv")

